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NEAR  FIELD  TURBULENT  WAKE  PREDICTIONS 


INTRODUCTION 

There  axe  several  experimental  and  numerical  studies  of  the  turbulent  wake  of 
a  self  propelled  body.  Schetz  and  Favin  (1977,  1979)  have  computed  the  near-body 
wake  using  an  elliptic,  axisymmetric  model  with  a  one-equation  turbulence  model. 
Comparisons  with  wind  tunnel  data  show  good  agreement  with  axial  velocity  but  an 
underprediction  of  the  swirl  velocities.  Swanson  and  Schetz  (1975)  made  numerical 
predictions  of  the  fax  wake  using  the  parabolic  axisymmetric  Navier- Stokes  equations 
with  a  one- equation  turbulence  model  at  locations  far  enough  downsteam  where  the 
swirl  velocities  could  be  neglected.  Swean  (1987)  has  simulated  the  turbulent  wake 
for  a  large  range  of  downstream  locations  including  axial  and  swirl  velocities  using 
the  three-dimensional  parabolic  Navier-Stokes  equations  and  two-equation  turbulence 
model  with  a  sophisticated  Reynolds  stress  closure.  In  general,  comparisons  with  the 
wind  tunnel  data  of  Mitra,  Neu,  and  Schetz  (1985,  1986)  were  good  especially  when  an 
account  was  made  of  the  axial  pressure  gradient  in  the  otherwise  parabolic  model. 

The  purpose  of  the  present  study  is  to  develop  a  partially  parabolic  model  of  the 
wake  of  a  self  propelled  body.  The  goal  is  to  eliminate  the  need  for  measured  down¬ 
stream  pressure  data  used  in  Swean’s  1987  study  while  maintaining  the  overall  good 
agreement  with  the  experimental  data.  The  following  sections  contain  a  description  of 
the  numerical  model.  A  description  of  the  experiments  and  Swean’s  (1987)  model  are 
available  in  his  report. 

NUMERICAL  MODEL 

Simplifying  Assumptions 

The  most  common  form  of  numerical  model  for  the  Navier-Stokes  equations  is 
elliptic.  Numerically,  this  makes  no  assumptions  about  the  importance  of  the  various 
terms.  It  does  require  the  storage  of  the  dependent  variables  over  the  entire  domain 
of  interest.  In  the  case  of  a  turbulent  wake,  this  includes  6  dependent  variables  (three 
fluid  velocities,  pressure,  turbulent  kinetic  energy  and  dissipation)  over  a  domain  that 
is  relatively  limited  in  the  transverse  coordinates  (perpendicular  to  the  principle  flow 
direction)  but  which  may  be  orders  of  magnitude  larger  in  the  axial  direction. 

To  overcome  the  problem  of  computer  storage  for  the  dependent  variables,  the 
parabolic  assumptions  axe  usually  made  whenever  axial  diffusion  terms  can  be  ignored 
and  the  axial  velocity  suffers  no  reversal  (see  Patankar  and  Spalding,  1972).  The 
solution  may  be  constructed  by  starting  from  the  location  of  the  initial  conditions  and 
marching  downstream  storing  the  values  of  the  dependent  variables  at  the  current  axial 
location  only.  The  advantages  of  this  procedure  are  greatly  reduced  computer  storage 
and  often  reduced  cpu  time. 
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Frequently  a  situation  occurs  where  the  parabolic  assumptions  axe  valid  except 
that  the  axial  pressure  gradient  is  important  (see  Pratap  and  Spalding,  1976).  The 
procedure  here  is  similar  to  the  purely  parabolic  case  except  that  the  pressure  must 
now  be  stored  over  the  entire  domain  and  multiple  sweeps  through  the  flow  field  are 
required  to  calculate  the  pressure  field.  The  advantage  of  this  procedure  is  again 
reduced  computer  storage  requirements,  hovever  execution  times  are  not  improved. 
This  procedure  is  referred  to  as  the  partially  parabolic  method  and  we  will  attempt  to 
adapt  it  to  model  the  wind  tunnel  data  of  Mitra  et  al.  (1985,  1986)  for  the  self-propelled 
body  with  a  simulated  free  surface. 

Governing  Equations 

The  equations  describing  the  partially  parabolic  form  of  the  Navier-Stokes  equa¬ 
tions  are  shown  in  Table  1.  The  distinguishing  feature  of  these  equations  are  that  all 
stress  gradients  in  the  axial  direction  except  for  pressure  are  neglected.  For  this  we 
require  that  1)  there  is  no  flow  reversal  in  the  axial  direction  and  2)  the  Reynolds 
number  is  high  enough  that  the  Reynolds  and  viscous  stresses  are  significant  only  in 
the  directions  normal  to  the  streamlines. 

Typically,  this  form  of  the  Navier-Stokes  equations  has  been  used  to  model  interior 
and  exterior  flows  for  which  pressure  influences  are  transmitted  from  a  downstream  lo¬ 
cation  to  an  upstream  location.  The  downstream  pressure  disturbance,  for  examples, 
may  be  produced  by  an  obstruction  in  a  duct  (confined  flow)  or  a  highly  curved  sur¬ 
face  (exterior  flow)  and  is  numerically  transmitted  upstream  through  successive  sweeps 
through  the  domain.  In  these  situations  the  initial  conditions  used  to  start  the  calcu¬ 
lation  are  the  same  as  for  parabolic  flows  i.e.  three  velocity  components  and  turbulent 
kinetic  energy  and  dissipation. 

The  sequence  of  steps  used  to  calculate  this  type  of  a  partially  parabolic  flow  field 
is  the  following: 

1)  The  three  dimensional  pressure  field  is  assigned  guessed  values,  usually  zero. 

2)  The  initial  values  of  the  u,  v,  and  w  components  of  velocity  are  advanced  to 
the  next  station  using  the  parabolic  procedure.  Existing  values  of  pressure  are  used 
to  evaluate  the  pressure  gradient  terms  in  the  momentum  equations.  Any  auxiliary 
equations  such  as  turbulence  or  temperature  are  advanced  at  this  time  also. 

3)  The  pressure  correction  equations  are  used  to  correct  the  values  of  pressure 
and  velocity  to  satisfy  continuity.  The  pressure  correction  equations  are  cast  in  three 
dimensional  form  but  upstream  and  downstream  values  of  the  correction  are  assumed 
zero.  The  pressure  is  stored  for  the  next  axial  sweep. 

4)  A  new  downstream  station  is  selected  and  the  momentum  and  continuity  equa¬ 
tions  are  solved  as  in  steps  2  and  3.  This  procedure  continues  until  the  last  downstream 
location  is  reached. 
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5)  A  new  axial  sweep  is  initiated  using  the  initial  conditions  and  the  three  di¬ 
mensional  pressure  field  calculated  from  the  previous  sweep.  Steps  2  through  4  are 
repeated  until  the  pressure  correction  or,  equivalently,  the  mass  residual  becomes  suf¬ 
ficiently  small. 

An  important  feature  of  the  procedure  is  that  the  pressure  disturbance  is  transmit¬ 
ted  upstream  by  the  calculation  of  the  axial  pressure  gradient  causing  an  axial  velocity 
gradient  and  mass  residual  which  is  then  eliminated  by  adjusting  the  upstream  value 
of  pressure  via  the  pressure  correction  equation.  This  does  not  leave  the  initial  pres¬ 
sure  conditions  intact.  For  the  situation  where  the  pressure  disturbance  has  its  source 
downstream  from  the  initial  plane,  this  presents  no  problem. 

Unfortunately,  this  is  not  the  case  in  the  current  situation.  The  pressure  distur¬ 
bance  is  caused  by  the  body  and  propeller  upstream  of  the  computational  domain.  The 
initial  pressure  distribution,  which  must  be  specified,  has  a  large  influence  on  the  flow 
field  close  to  the  body  and  propeller. 

When  the  stand.ird  method  of  Pratap  and  Spalding  (1976)  is  used,  the  solution 
converges  relatively  quickly  with  small  mass  and  momentum  residuals.  However,  the 
procedure  has  no  mechanism  to  incorporate  the  effects  of  an  initial  pressure  distribu¬ 
tion  in  the  converged  solution.  The  result  is  that  the  initial  pressure  distribution  is 
annihilated  by  the  pressure  correction  equation  leaving  a  pressure  distribution  which 
does  no  work  on  the  flow  field.  The  velocity  fields  tend  to  relax  monotonically  toward 
free  stream  values  unlike  the  trends  apparent  in  the  wind  tunnel  measurements. 

Mathematically  it  is  unusual  to  find  both  the  normal  velocity  and  the  pressure  » 

specified  at  a  boundary  as  it  is  here.  The  problem  becomes  over-specified  at  that  lo¬ 
cation  and  there  are  not  enough  degrees  of  freedom  to  satisfy  both  the  Navier- Stokes 
equations  and  continuity  exactly.  The  standard  method  of  Pratap  and  Spalding  ( 1976) 
formulates  the  problem  with  non-zero  initial  conditions  on  velocity  only.  The  comple¬ 
mentary  formulation  using  the  initial  o-essure  distribution  instead  of  the  axial  velocity  • 

was  also  computed  with  a  similar  lack  of  success.  Instead  of  annihilating  the  initial 
pressure  as  before,  it  is  the  initial  axial  velocity  which  becomes  distorted. 

There  is  a  method  in  the  literature  which  naturally  incorporates  the  Dirichlet  type 
of  boundary  condition  on  the  pressure  equation  at  the  initial  plane.  Where  the  original 
method  calculates  the  elliptic  pressure  field  in  a  marching  procedure  starting  at  the  *- 

initial  plane  and  marching  downstream,  a  method  similar  to  Chilukuri  and  Pletcher 
(1980)  can  be  used  which  solves  the  three-dimensional  pressure  equation  and  leaves  the  j 

initial  pressure  conditions  intact.  The  basic  assumption  here  is  that  either  the  initial  , 

pressure  gradient  is  zero  (which  again  eliminates  the  influence  of  the  initial  pressure) 
or,  the  initial  pressure  and  velocity  conditions  satisfy  the  Navier-Stokes  equations.  * 

the  continuity  equation,  and  the  modeled  form  of  the  Reynolds  stresses  (which  is  the 
assumption  used  here).  Mathematically,  the  difference  between  these  two  methods 
is  that  the  Pratap  method  solves  an  approximate  equation  which  is  used  to  correct 
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a  previous  estimate  of  pressure  where  the  method  of  Chilukuri  calulates  the  actual 
elliptical  pressure  field  based  on  values  of  the  velocities  using  the  full  Navier-Stokes 
equations. 

The  goal  here  is  to  find  a  numerical  procedure  using  either  form  of  the  pres¬ 
sure/continuity  equation  which  compromises  between  satisfying  the  Navier-Stokes 
equations  and  continuity  while  incorporating  both  the  initial  axial  velocity  distribu¬ 
tion  and  initial  pressure.  This  compromise  will  only  be  noticable  near  the  initial  plane 
and  will  depend  on  the  quality  of  the  initial  conditions. 

Method  Development 

As  mentioned  before,  the  first  method  used  was  the  procedure  of  Pratap  and 
Spalding  (1972).  Because  this  method  is  unable  to  use  initial  pressure  conditions  the 
solution  converges  to  a  mathematically  satisfying  but  physically  unrealistic  flow  field. 
In  order  to  correct  this  deficiency  in  the  method,  a  modification  was  used  to  keep  the 
initial  pressure  field  intact.  At  the  beginning  of  each  sweep,  the  axial  velocity  field 
was  advanced  using  the  pressure  gradient  between  the  initial  pressure  plane  and  the 
first  calculated  pressure  plane  downstream.  The  transverse  velocity  fields  are  advanced 
to  the  next  downstream  station.  The  pressure  correction  equation  is  used  to  correct 
the  transverse  velocity  components  but  not  the  initial  pressure.  At  the  second  and 
all  subsequent  stations  downstream,  the  velocity  fields  are  advanced  and  the  pressure 
correction  equation  is  used  to  correct  both  the  pressure  and  transverse  velocties. 

This  represents  a  simplistic  way  to  overcome  the  initial  pressure  field  problem. 
Unfortunately,  the  solution  diverges  at  the  initial  plane  location.  The  connection  be¬ 
tween  the  pressure  at  the  initial  plane  and  the  first  station  downstream  is  too  weak 
and  a  discontinuity  in  the  steamwise  variation  of  the  axial  velocity  and  pressure  fields 
results. 

To  strengthen  the  steamwise  connections  in  the  pres°*Te  field,  several  modifications 
were  used  but  without  much  success.  A  three  dimensional  elliptic  version  of  the  pressure 
correction  equation  was  solved  at  all  axial  stations  including  the  initial  plane.  Here 
again,  the  initial  correction  was  used  only  to  correct  the  velocity  field  leaving  the  initial 
pressure  condition  intact.  This  calculation  also  diverged.  Instead  of  simply  forming 
a  discontinuity  in  the  axial  velocity  field  as  the  previous  attempt  did.  this  calculation 
eventually  became  unstable  again  at  the  location  of  the  initial  plane. 

At  a  very  minor  increase  in  core  requirements  over  the  three  dimensional  pressure 
correction  algorithm,  the  method  of  Chilukuri  et  al.  (1980)  can  be  used.  This  method 
does  not  use  the  pressure  correction  equation.  Instead,  the  discretized  forms  of  the 
Navier-Stokes  equations  are  used  in  the  continuity  equation.  After  some  manipulation, 
the  equations  can  be  separated  into  a  standard  form  for  the  full  pressure  equations  with 
source  terms  containing  the  velocities  from  the  previous  sweep.  Although  this  source 
term  appears  to  be  difficult  to  compute  accurately,  it  actually  is  already  available  from 


4 


the  previous  calculations  of  the  velocity  equations.  In  this  way,  there  is  little  increase 
in  computational  effort  or  memory  compared  to  the  3-D  pressure  correction  method. 
The  overall  structure  of  the  program  is  affected  because  the  two  dimensional  pressure 
correction  equation  must  be  retained  both  for  the  initial  sweep  of  the  entire  domain  to 
generate  the  initial  guess  of  the  downstream  pressure  field  and  to  maintain  continuity 
in  the  regions  where  the  axial  pressure  gradient  is  ignored. 

There  are  two  potential  advantages  using  this  method.  The  first  one,  which  is 
often  found  in  the  literature,  is  that  the  convergence  rate  is  improved  because  the  full 
Navier-Stokes  equations  are  used  to  satisfy  continuity  rather  than  the  approximate  form 
of  the  pressure  correction  method.  More  important  for  the  current  problem  is  that  the 
boundary  conditions  for  the  elliptic  equation  are  in  the  form  of  pressure  rather  than 
a  pressure  correction  which  must  be  surpressed  at  the  initial  plane.  The  coefficients 
linking  the  pressure  field  to  the  initial  pressure  conditions  are  easily  calculated.  It  does 
not,  however,  alter  the  over-specified  condition  at  the  initial  plane. 

In  practice  the  faster  pressure  convergence  rate  expected  of  this  method  was  ap¬ 
parent  over  most  of  the  flow  field.  It  did  not  occur  at  the  first  axial  station  and,  in 
fact,  the  calculation  slowly  diverged  and  eventually  became  unstable  at  that  location. 

The  last  two  methods  which  have  been  described  involve  three  dimensional  calcu¬ 
lation  of  the  pressure  field  and,  as  a  result,  require  considerably  more  computer  storage. 
The  memory  limitations  of  the  LCP2  Vax  11-780  threaten  to  make  the  axial  grid  spac¬ 
ing  unacceptably  coarse.  An  examination  of  the  parabolic  calculations  showed  that  the 
effect  of  the  axial  pressure  gradient  became  negligible  at  about  x/D  =  6.  Thus  the  pro¬ 
gram  was  modified  so  that  the  partially  parabolic  region  was  confined  to  0  <  x/D  <  6 
and  a  parabolic  calculation  was  used  to  calculate  the  rest  of  the  flow  field  6  <  x/D  < 
28.  An  additional  advantage  to  this  procedure  is  that  the  multiple  sweeps  are  restricted 
to  the  reduced  domain  and  the  full  domain  is  calculated  only  on  the  last  sweep,  greatly 
reducing  cpu  requirements. 

Modified  Flow  Models 

All  attempts  to  make  the  flow  field  conform  to  all  parts  of  the  model  within  the 
first  step  downstream  creates  a  highly  distorted  pressure  field  which,  in  turn,  causes 
the  axial  velocity  gradient  to  become  distorted.  The  basic  method  which  appears  to 
avoid  this  instability  is  to  separate  the  calculation  of  the  axial  pressure  gradient  from 
the  transverse  pressure  gradient  in  some  way.  Two  approaches  were  used  and  their 
advantages  and  disadvantages  will  be  discussed. 

Pressure  Correction  Method 

There  aLre  several  advantages  to  using  the  pressure  correction  method  in  modified 
form  for  this  problem.  The  most  obvious  is  that  computer  storage  requirements  are 
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small.  The  second  is  that  all  disturbances  in  this  type  of  flow  situation  are  propogat- 
ing  downstream.  This  eliminates  one  of  the  primary  advantages  of  three  dimensional 
pressure  calculation  methods  which  is  that  downstream  pressure  disturbances  are  felt 
throughout  the  domain  in  one  axial  sweep.  The  final  advantage  is  that  the  equation 
is  directly  formulated  to  annihilate  the  mass  residual.  It  uses  an  approximate  form 
of  the  Navier-Stokes  equations  in  the  continuity  equation  and  has  the  property  that 
the  calculated  pressure  field  will  satisfy  the  Navier-Stokes  equations  only  in  the  limit 
of  zero  mass  residual.  The  approximate  nature  of  the  initial  conditions  should  be  less 
important  using  this  method. 

The  marching  form  of  the  pressure  correction  equation  is  used  here.  During  each 
sweep,  the  pressure  is  stored  at  each  axial  downstream  location  but  the  pressure  gradi¬ 
ent  in  the  axial  velocity  equation  is  used  from  the  previous  sweep.  Notice  that  by  using 
the  marching  version  of  the  pressure  correction  equation  we  maintain  the  advantage  of 
reduced  computer  storage  requirements  compared  to  full  three  dimensional  pressure  or 
pressure  correction  methods. 

In  this  way,  a  purely  parabolic  calculation  is  used  to  provide  the  initial  axial 
pressure  distribution.  All  subsequent  sweeps  will  use  the  calculated  pressure  gradient. 
The  difference  between  this  formulation  and  the  standard  pressure  correction  method 
is  that  the  initial  guess  for  the  transverse  pressure  field  at  each  axial  station  is  advanced 
from  the  immediate  upstream  location.  The  difference  between  the  marching  scheme 
and  the  three-dimensional  methods  is  that  axial  effects  in  the  marching  scheme  are 
included  in  the  pressure  equation  only  through  the  neighboring  pressure  correction 
coefficients.  The  values  of  actual  upstream  and  downstream  pressure  corrections  are 
assumed  to  be  zero.  Only  the  transverse  velocities  are  adjusted  using  the  pressure 
corrections.  If  the  axial  velocity  is  also  adjusted,  the  effect  is  the  same  as  modifying 
the  axial  pressure  gradient  which  effectively  removes  the  upstream  pressure  conditions 
and  the  resulting  pressure  gradient  evenutally  relaxes  to  a  version  which  again  does  no 
work  on  the  axial  velocities. 

By  advancing  the  upsteam  transverse  pressure  as  an  initial  estimate  and  correcting 
only  the  transverse  velocities,  the  effect  of  the  initial  pressure  conditions  is  transmit¬ 
ted  downstream  because  each  new  transverse  pressure  field  is  exactly  the  sum  of  the 
upstream  pressure  plus  the  correction  needed  for  the  actual  mass  residual  caused  by 
advancing  one  axial  location.  At  each  axial  location,  only  one  iteration  of  the  pressure 
equation  is  needed.  Because  only  the  transverse  velocities  are  corrected  during  the 
sweep,  the  mass  residual  cannot  be  eliminated  during  the  current  sweep,  rather,  sev¬ 
eral  sweeps  are  needed  so  that  the  axial  velocity  may  respond  to  the  new  pressure  field. 
Several  implicit  iterations  at  the  same  axial  location  would  produce  a  slightly  smaller 
mass  residual  but  subsequent  pressure  corrections  would  not  be  accurate.  The  first 
iteration  would  produce  a  pressure  and  velocity  correction  based  on  the  actual  mass 
residual.  Further  iterations  at  the  same  location  would  produce  pressure  corrections 


based  on  the  axial  part  of  the  mass  residual  which  cannot  be  immediately  eliminated. 
A  distorted  pressure  field  would  result.  During  the  last  axial  sweep,  however,  multiple 
implicit  iterations  should  be  used  as  the  pressures  calculated  then  will  not  be  used  for 
the  axial  pressure  gradient  and  the  greatest  possible  reduction  of  the  mass  residual  is 
needed  to  accurately  predict  swirl  velocities. 

Similarly,  the  mass  residual  could  be  eliminated  during  a  single  sweep  by  writing 
the  pressure  correction  equation  only  for  the  transverse  coordinates  and  eliminating  the 
upstream  and  downstream  coefficients.  This  forces  the  transverse  velocities  to  be  solely 
responsible  for  satisfying  continuity.  The  result  is  that  at  the  first  station  downstream 
of  the  initial  condition  plane  the  maximum  swirl  velocity  was  almost  doubled  and  the 
calculated  transverse  pressure  was  no  longer  related  to  the  axial  variation  of  pressure 
and  could  not  be  used  for  this  on  the  next  sweep. 

The  final  result  of  this  calculation  was  a  flow  field  which  converged  within  15  axial 
sweeps.  Near  the  initial  condition  plane,  a  finite  mass  residual  remains  but  downstream 
it  quickly  dies  out.  This  finite  residual  appears  to  be  caused  by  the  approximate  nature 
of  the  initial  conditions.  Comparisons  with  the  experimental  data  will  be  presented  in 
a  later  section. 

Total  Pressure  Formulation 

This  method  is  actually  a  small  modification  of  the  method  of  Chilukuri  ( 1980 ) 
and  uses  the  same  general  principles  as  the  previous  method.  Here,  the  axial  sweep 
of  the  velocities  and  transverse  pressure  is  exactly  the  same.  On  the  first  sweep,  the 
calculation  is  again  purely  parabolic  with  no  axial  pressure  gradient  and  the  resulting 
three-dimensional  pressure  field  is  stored  as  an  initial  guess.  On  subsequent  sweeps,  the 
pressure  correction  equation  is  calculated  to  satisfy  the  transverse  continuity  condition 
in  the  same  way  as  the  previous  method.  The  pressure  field  calculated  bv  the  pressure 
correction  method  is  not  stored.  Instead,  the  three  dimensional  coefficients  for  the  full 
pressure  equation  are  calculated  and,  at  the  completion  of  each  axial  sweep,  the  full 
pressure  field  is  calculated. 

Here  since  we  are  not  using  the  values  of  the  calculated  transverse  pressure  field 
for  the  axial  pressure  gradient,  multiple  implicit  iterations  at  each  axial  location  are 
useful  to  provide  the  best  estimate  of  the  converged,  mass  conserving  velocity  fields 
using  the  current  estimate  of  the  axial  pressure  variation.  This  will  help  provide  the 
best  estimates  of  the  coefficients  for  the  pressure  equation. 

This  method  has  some  of  the  same  disadvantages  of  the  original  version  which  is 
increased  computer  storage  over  the  pressure  correction  method.  In  addition,  because 
a  hybrid  procedure  is  used,  i.e.  the  pressure  correction  method  is  solely  used  for 
continuity  and  velocity  correction  and  the  total  pressure  method  is  used  to  calculate 
the  axial  pressure  gradient,  the  calculation  converges  much  more  slowly  requiring  about 
100  iterations.  In  addition,  this  method  uses  multiple  implicit  iterations  on  all  but  the 
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first  axial  sweep  while  the  pressure  correction  method  uses  multiple  iterations  only 
during  the  last  sweep.  The  final  result  is  that  the  pressure  correction  method  required 
approximately  1.75  hours  of  cpu  on  the  LCP2  vax  while  the  total  pressure  method 
required  approximately  25  hours.  How  much  of  this  difference  was  due  to  increased 
calculations  and  how  much  was  caused  by  the  increased  need  for  memory  paging  due 
to  the  much  larger  memory  requirements  is  unknown.  Of  course  either  method  would 
require  negligable  cpu  times  on  the  NRL  Cray. 

The  advantage  of  this  method  is  that  more  of  the  physics  of  the  problem  are 
included  in  this  calculation  of  the  axial  pressure  variation  because  the  form  of  the 
equation  is  nominally  exact.  As  we  shall  see,  the  differences  between  the  predictions  of 
the  two  methods  tire  relatively  small. 

COMPARISONS  WITH  EXPERIMENTS 

The  experiments  used  to  compare  with  the  numerical  predictions  have  been  ex¬ 
tensively  reported  in  Mitra,  Neu,  and  Schetz  (1985,  1986).  The  experiments  were 
essentially  of  double  axisymmetric  bodies  connected  by  a  thin  strut  in  a  wind  tunnel. 
They  were  used  to  simulate  a  SWATH  (Small  Waterplane  Area  Twin  Hull)  type  ship 
and  the  plane  of  symmetry  represents  the  free  surface.  No  attempt  was  made  to  sim¬ 
ulate  the  second  half  of  this  type  of  ship  as  the  hulls  are  separated  far  enough  in  an 
actual  ship  that  the  merging  of  the  wakes  would  not  occur  in  the  domain  of  interest. 
Several  sets  of  experiments  were  run.  In  the  situation  we  will  be  concerned  with  here, 
the  free  stream  velocity  was  45.72  m/s,  both  the  body  and  its  mirror  image  were  self 
propelled  with  a  15.24  cm  propeller  with  thrust  balanced  against  drag  using  a  strain 
guage,  and  transverse  plane  data  was  taken  at  3.17  cm  (x/D  =  .208),  15.24  cm.  (x/D 
=  1),  and  427  cm.  (x/D  =  28). 

In  addition,  a  detailed  comparison  between  these  measurements  and  a  parabolic 
finite  element  code  is  reported  in  Swean  (1987).  In  Swean's  report  comparisons  were 
made  with  the  mean  velocities,  turbulent  kinetic  energy,  and  Reynolds  stresses.  In 
order  to  accurately  model  the  self-propelled  body,  Swean  used  the  measured  pressure 
profiles  at  the  three  axial  positions  to  estimate  the  axial  pressure  gradient.  In  general, 
the  agreement  was  very  good.  The  purpose  of  the  current  study  was  to  present  a 
calculations!  model  which  would  make  predictions  based  solely  on  initial  plane  data. 

Here  we  will  compare  the  predictions  of  the  present  models  with  the  experimental 
data  and  the  previous  parabolic  predictions.  Only  those  features  of  the  data  which 
axe  influenced  by  the  partially  parabolic  features  of  the  code  will  be  presented.  The 
parabolic  finite  element  code  used  by  Swean  solved  the  same  set  of  equations  describing 
momentum  transport  and  turbulent  kinetic  energy  and  dissipation  as  the  present  finite 
volume  code.  Continuity  was  satisfied  in  both  using  a  Poisson  type  equation  but  the 
current  formulation  uses  a  predictor-corrector  type  formulation  while  the  finite  element 
code  uses  a  penalty  method.  To  differentiate  between  the  two  methods  and  a  purely 
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parabolic  (P)  calculation,  we  will  refer  to  Swean’s  method  as  a  quasi-partiallv  parabolic 
model  (QPP)  and  the  present  method  as  modified  partially  parabolic  (MPP' 

Even  though  experimental  data  is  available  to  initiate  the  calculation,  this  data 
covers  only  part  of  the  domain  (see  Swean,  1987  for  details).  The  data  was  interpolated, 
extrapolated, and  reflected  to  cover  the  missing  areas.  This  was  made  especially  difficult 
because  of  the  absence  of  true  symmetry  due  to  the  rotation  of  the  fluid  caused  by  the 
propeller.  This  process  contributes  to  the  approximate  nature  of  the  initial  conditions. 

Here,  we  use  the  initial  conditions  as  formulated  by  Swean  including  the  grid,  the 
velocity  fields,  the  turbulent  kinetic  energy,  and  the  turbulent  dissipation.  The  only 
modification  made  was  to  fit  the  data  onto  a  staggered  grid. 

Comparisons  between  the  MPP  and  QPP  predicted  and  measured  values  of  the 
axial  minimum  and  maximum  velocity  are  presented  in  figures  1  a  and  b.  The  most 
striking  feature  of  these  plots  is  the  peak  of  the  maximum  velocity  which  occurs  down¬ 
stream  of  the  initial  plane.  This  peak  is  caused  by  the  axial  pressure  gradient  and  shows 
the  need  for  the  partially  parabolic  assumptions.  It  is  interesting  that  the  predicted 
axial  location  of  this  peak  and  the  location  of  the  measured  values  coincides.  No  special 
effort  was  used  in  the  calculation  to  insure  this  feature.  A  strictly  parabolic  calculation 
with  no  streamwise  pressure  gradient  does  not  predict  a  downstream  maximum  for  the 
axial  velocity. 

The  QPP  prediction  of  this  data  set  also  showed  that  the  measured  and  predicted 
peak  occured  at  the  same  axial  location  but  this  was  due  to  the  use  of  the  measured  axial 
pressure  gradient  which  places  the  peak  at  that  location.  The  value  of  the  measured 
peak  axial  velocity  is  greater  than  the  predicted  value  in  both  the  QPP  and  MPP 
predictions.  In  the  MPP  calculations,  the  underprediction  appears  to  be  caused  by 
the  initial  adjustment  region  directly  downstream  of  the  initial  plane.  During  these 
initial  calculations,  mass  continuity  is  not  perfectly  satisfied  and  some  of  the  effect  of 
the  pressure  gradient  appears  to  be  lost. 

The  QPP  calulation  requires  that  continuity  be  satisfied  to  within  a  chosen  toler¬ 
ance  before  advancing  to  the  next  axial  station  yet  the  peak  value  of  the  axial  \-elocity 
is  less  than  either  the  measured  values  or  the  current  predictions.  There  are  two  pos¬ 
sible  explanations  for  this.  The  simplest  reason  is  that  it  uses  a  value  of  the  pressure 
gradient  which  is  an  average  between  the  initial  plane  and  x/D  =  1  downstream  and  | 

which  has  been  extensively  extended  to  cover  the  entire  transverse  plane.  The  other 
possible  cause  may  be  due  to  the  method  in  which  overall  continuity  is  satisfied.  For 
the  present  MPP  models,  during  any  sweep  through  the  axial  domain,  the  axial  pres¬ 
sure  gradient  used  was  calculated  from  the  previous  sweep.  The  traditional  method  of 
satisfying  continuity  at  each  axial  location  is  to  correct  all  three  components  of  veloc- 
ity.  It  is  easy  to  show  that  by  correcting  the  axial  velocity  in  the  MPP  method,  the 
effect  is  to  change  (usually  reducing)  the  axial  pressure  gradient.  Assuming  for  the 
moment  that  the  finite  element  QPP  code  is  similar  in  this  respect,  the  downstream 
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axial  velocity  peak  would  have  been  reduced  by  a  strict  continuity  requirement  at  the 
first  several  stations.  This  could  possibly  explain  the  smaller  peak  found  in  the  QPP 
model  compared  to  the  MPP  model. 

The  other  profile  shown  is  the  minimum  axial  velocity  either  at  the  body  center 
or  behind  the  strut  at  the  simulated  free  surface.  Initially  this  minimum  occurs  at  the 
center  of  the  body  but,  at  some  distance  past  the  x/D=l  station,  the  velocity  deficit 
behind  the  strut  at  the  symmetry  line  becomes  the  minimum.  The  use  of  the  axial 
pressure  gradient  in  both  the  QPP  and  MPP  calculations  causes  the  velocity  deficit 
in  both  locations  to  relax  back  to  free  stream  values  much  more  quickly  than  a  purely 
parabolic  prediction.  In  fact,  the  current  MPP  as  well  as  the  QPP  calculation  appears 
to  have  slightly  overpredicted  the  relaxation  rate  compared  to  the  experiments. 

Qualitatively,  both  the  QPP  and  MPP  calculations  show  the  same  behavior.  Quan¬ 
titatively,  the  velocity  defect  behind  the  body  center  for  the  MPP  method  continues 
to  relax  back  to  the  free  stream  value  faster  than  the  QPP  method.  This  changes  the 
location  where  the  minimum  axial  velocity  moves  from  the  body  center  to  the  intersec¬ 
tion  of  the  strut  and  symmetry  line.  For  the  MPP  method,  this  occurs  at  about  x/D 
=  3  and  for  the  QPP  method  at  about  x/D  =  7.  The  experimental  measurements  only 
shows  that  it  occurs  sometime  between  the  x/D  =  1  and  x/D  =  28  stations. 

Figures  2,  3,  and  4  shows  the  development  of  the  swirl  velocity  compared  to  the 
experimental  values.  In  the  plot  of  values  at  x/D  =  1  the  predictions  agree  well  with  the 
measured  values.  The  maximum  value  of  swirl  is  slightly  less  than  either  the  measured 
or  the  QPP  model.  The  reason  for  this  is  that  this  value  is  very  sensitive  to  the  mass 
continuity  calculations  upstream.  Since  continuity  is  not  used  as  a  criterion  for  the 
local  convergence  during  each  sweep,  a  set  number  of  iterations  (from  3  to  10)  were 
performed  at  each  axial  location  before  advancing  downstream.  The  swirl  predictions 
improved  as  the  number  of  iterations  increased  but  the  overall  execution  time  quickly 
became  excessive.  The  calculations  using  20  iterations  in  the  last  sweep  of  the  pressure 
correction  form  of  the  MPP  model  showed  that  this  minor  difference  (seen  here  for  10 
iterations)  was  due  to  this  sensitivity. 

A  more  interesting  difference  is  seen  in  the  comparison  at  x/D  =  28.  Again  the 
value  for  the  maximum  swirl  velocity  calculated  by  the  MPP  method  is  slightly  low. 
The  center  of  the  swirl  has  also  moved  slightly  to  the  right  unlike  the  QPP  method. 
The  reason  for  this  illustrates  the  differences  in  the  computational  methods.  The  finite 
volume  calculation  uses  a  hybrid  upwind/ central  difference  method.  The  finite  element 
method  uses  something  similar  to  artificial  dissipation  to  maintain  numerical  stability 
and  the  effect  of  this  was  minimized  for  this  calculation  by  Swean.  Thus  the  finite 
volume  calculation  has  considerably  more  numerical  diffusion  and  the  effect  of  the 
mirror  image  body  and  propeller  was  felt  more  quickly.  The  effect  of  this  image  is  an 
induced  translation  of  the  propeller  vortex  in  the  direction  of  the  tangential  velocity 
at  the  symmetry  line  similar  to  the  self  induced  translation  of  a  vortex  pair  (see  for 


example  Sarpkaya,  1982).  While  the  effect  is  real,  it  is  obviously  overestimated  here 
due  to  the  excessive  numerical  diffusion.  A  simple  solution  for  this  problem  is  the  use 
of  third  order  upwind  techniques  such  as  the  QUICK  scheme  (Freitas  et  al.  1985).  At 
this  axial  location  some  countersign  vorticity  begins  to  appear  where  the  swirl  velocities 
interact  with  the  right  boundary.  This  is  again  caused  by  the  translation  of  the  propeller 
vortex.  The  predictions  are  good  nevertheless. 

A  comparison  of  the  axial  decay  rate  of  turbulent  kinetic  energy  at  two  locations  for 
both  the  MPP  and  QPP  calculations  is  shown  in  figures  5  a  and  b.  The  top  line  shows 
k  for  the  location  directly  behind  the  body  at  the  propeller  center.  At  the  initial  plane, 
the  body  center  has  a  relatively  low  level  of  turbulent  energy  surrounded  by  higher 
levels.  The  QPP  method  predicts  that  at  the  x/D  =  1  axial  station,  the  surrounding 
turbulent  energy  was  expected  to  diffuse  inward,  eliminating  the  depression  in  k  and 
raising  the  value  above  that  at  the  initial  plane.  Measurements  seem  to  indicate  much 
less  diffusion  was  present  leaving  the  depressed  center  area  intact  with  a  value  less  than 
the  initial  value  at  x/D  =  1.  As  shown  here,  the  MPP  predictions  are  that  the  value  at 
the  body  center  does  drop  going  downsteam.  This  agreement  with  the  measurements 
may  be  more  fortuitous  than  accurate  however.  In  figures  6  a-d  the  radial  distribution 
of  k  at  x/D  =  1  is  shown.  Although  the  center  value  is  quite  reasonable,  the  value 
of  the  surrounding  peak  is  depressed.  The  figure  showing  the  QPP  results  at  that 
location  show  that  although  the  center  depression  is  eliminated,  the  estimates  of  the 
surrounding  peak  are  quite  good.  It  appears  that  in  the  current  model  the  generation 
terms  in  the  turbulent  kinetic  energy  equations  are  too  crudely  estimated. 

The  original  programming  was  designed  to  produce  readily  recognizable  code,  thus' 
the  various  Reynolds  stress  components  which  went  into  the  generation  terms  were 
left  intact.  Derivatives  of  the  Reynolds  stresses  were  made  using  the  discrete  volume 
values  which  were  in  turn,  found  from  the  surrounding  values  of  velocity,  k  and  e. 
The  result  is  that  a  point  value  for  the  generation  term  is  calculated  from  distant 
values  of  the  other  primitive  variables  and  the  smea.  hig  causes  an  underestimate.  In 
the  finite  element  code,  all  derivatives  are  evaluated  at  the  point  they  are  used  with 
second  derivatives  reduced  to  products  of  first  derivatives.  The  result  is  that  much  less 
diffusion  is  expected. 

A  future  improvement  on  this  code  would  be  to  reduce  the  source  terms  in  the 
turbulence  equations  to  primitive  variables  (a,  v,  w ,  k,  and  e).  This  would  require  some 
fairly  tedious  manipulations  and  the  coding  would  be  much  less  transparent  but  nu¬ 
merical  diffusion  should  be  reduced  and  the  peaks  in  turbulent  kinetic  energy  should  be 
higher.  Whether  the  center  depression  in  turbulent  kinetic  energy  would  be  maintained 
is  not  clear  and  the  numerical  diffusion  would  still  be  greater  than  that  found  in  the 
finite  element  calculations. 

The  experimental  data  has  a  second,  much  higher  peak  value  of  turbulent  energy 
at  about  R/D  =  .4  as  shown  in  figure  6.  This  peak  actually  increases  by  about  50%  by 


x/D  =  1  from  the  initial  plane.  All  of  the  models  predicted  approximately  equivalent 
results  for  this  feature.  Instead  of  increasing,  the  peaks  decayed  by  about  50%.  At  this 
location,  velocity  gradients  and  predicted  turbulent  production  are  small.  It  appears 
the  the  measured  increase  in  turbulent  energy  is  actually  caused  by  velocity  fluctuations 
due  to  the  passing  of  the  propeller  tip.  An  analysis  of  the  turbulent  spectrum  by  Mitra 
et.  al.  (1986)  shows  a  large  peak  at  the  primary  and  harmonic  frequencies  associated 
with  propeller  rotation. 

Comparisons  Between  Current  Methods 

Thus  far,  all  results  shown  have  been  for  the  total  pressure  formulation  version 
of  this  model.  As  was  described  earlier,  this  version  uses  the  complete  Navier- Stokes 
equtions  in  the  continuity  equation  to  calculate  pressure.  It  was  used  primarily  because 
it  contains  more  of  the  physics  in  complete  form  than  does  the  method  using  the 
pressure  correction.  There  are  two  reasons  for  preferring  the  pressure  correction  version 
if  the  accuracy  is  comparable.  First,  the  execution  time  is  reduced  by  a  factor  of  16 
and,  second,  storage  is  reduced  by  a  factor  of  8.  Of  course,  the  magnitude  of  these 
reductions  are  specific  to  this  problem  but  they  appear  to  be  conservative  estimates. 

Some  general  observations  concerning  the  predictions  produced  by  both  versions 
will  illustrate  the  differences.  First,  they  differ  primarily  in  the  way  they  calculate  the 
axial  pressure  distribution.  The  total  pressure  version  can  be  expected  to  produce  a 
better  estimate  of  the  axial  pressure  distribution  at  the  expense  of  small  errors  in  the 
satisfaction  of  continuity  near  the  initial  plane.  This  is  because  the  axial  pressure  dis¬ 
tribution  is  calculated  using  a  global,  three  dimensional  calculation  but  the  equations 
used  axe  different  from  the  ones  used  to  calculate  transverse  pressure  to  satisfy  conti¬ 
nuity.  The  pressure  correction  method  is  expected  to  satisfy  the  continuity  equation 
better  but  produce  a  cruder  estimate  of  the  axial  pressure  gradient.  This  is  because 
axial  pressure  is  calculated  using  an  approximate  equation  but  the  axial  and  trans¬ 
verse  pressure  calculations  are  consistant  (if  separate).  Finally,  if  continuity  is  satisfied 
exactly,  both  methods  should  give  equivalent  results. 

Figure  7  contains  a  comparison  between  the  decay  of  the  maximum  swirl  velocity 
for  the  total  pressure  version,  the  pressure  correction  version,  and  the  measurements. 
The  pressure  correction  version  is  better  at  satisfying  the  transverse  part  of  the  problem. 
At  x/D  =  1  the  maximum  swirl  is  equivalent  to  the  QPP  finite  element  code  and  closer 
to  the  experimental  measurements.  In  addition,  the  maximum  mass  residual  (which 
naturally  occurs  just  downstream  of  the  initial  conditions)  is  much  less  for  the  pressure 
correction  method  than  the  total  pressure  version. 

The  total  pressure  version  appears  to  be  better  at  satisfying  the  axial  part  of  the 
problem.  Comparisons  between  the  maximum  and  minimum  axial  velocity  as  shown 
in  figure  8  are,  in  general,  better  using  this  version.  In  addition,  mass  continuity  is 
satisfied  much  better  using  this  method  after  a  few  downstream  steps.  This  is  expected 


once  the  problems  associated  with  the  initial  conditions  die  out.  After  about  x/D  =  5, 
both  methods  give  equivalent  results.  This  is  expected  since  the  mass  residual  of  both 
quickly  becomes  small  downstream  of  the  initial  plane. 

Comparisons  between  predicted  and  measured  pressure  at  the  downstream  loca¬ 
tions  show  only  general  similarities.  Figure  9  shows  the  variation  of  the  minimum  and 
maximum  pressure  as  a  function  of  axial  position.  The  trends  are  the  same  but  the 
maximum  pressure  difference  is  larger  for  the  pressure  correction  method  and  appears 
to  be  closer  to  the  measured  values. 

Figures  10  a  -  c  shows  that  the  calculated  transverse  pressure  distribution  had 
only  a  general  resemblence  to  the  measured  distribution  at  x/D  =  1.  Figure  10a  shows 
the  interpolated  experimental  pressure  data  at  x/D  =  1  and  figures  10  b-c  show  the 
distributions  computed  by  each  of  the  two  methods  described  above.  With  reference  to 
figure  10a,  however,  it  should  be  noted  that  the  experiments  were  principally  confined  to 
the  2nd  quadrant  (upper  left)  with  the  only  additional  data  existing  along  the  9  =  225° 
and  270°  radials  (as  measured  counter-clockwise  from  the  right  horizontal  axis).  Since 
the  only  natural  plane  of  symmetry  for  the  wind-tunnel  configuration  is  z/D  =  1,  it  was 
necessary  to  interpolate  very  sparse  data  to  initialize  the  3rd  quadrant  after  which  the 
data  were  reflected  to  quadrants  4  and  1,  either  symmetrically  or  anti-symmetrically  as 
appropriate.  As  a  consequence  any  data  structure  shown  can  be  misleading  especially  in 
quadrants  3  and  4.  The  primary  development  in  the  calculations  here  is  the  transition 
of  the  pressure  at  the  center  of  the  body  from  a  large  positive  pressure  to  the  minimum 
value  of  negative  pressure.  Surrounding  the  body  center  is  a  region  of  high  pressure. 
Otherwise,  the  experimental  values  show  much  more  detailed  structure.  Whether  or  not 
some  of  this  structure  is  the  result  of  the  extrapolation  and  interpolation  routines  used 
to  map  the  data  over  the  entire  domain  is  unknown.  Here  again  there  are  differences 
between  the  two  calculation  methods.  The  complete  pressure  formulation  captured 
the  tighter  structure  better  while  the  pressure  correction  method  estimated  the  total 
pressure  difference  better. 

In  figures  11  a  and  b,  the  calculated  pressure  distribution  at  x/D  =  28  definitely 
shows  the  influence  of  the  image  body  and  the  drag  wake  of  the  strut  on  the  calculation. 
Interestingly,  this  is  also  clearly  seen  in  the  experimental  data.  At  this  axial  location 
there  was  a  strong,  rectangular  depression  in  the  experimental  pressures  which  seems 
to  be  either  an  artifact  of  the  measurement  technique  or  a  pressure  drop  due  to  the 
wind  tunnel  walls  and  this  was  removed.  As  was  described  earlier,  the  effect  of  the 
image  body  was  not  evident  in  the  swirl  velocity  data  at  this  axial  location. 

CONCLUSIONS 

The  near  field  turbulent  wake  of  a  self  propelled  body  presents  several  interesting 
challenges  for  numerical  prediction  schemes.  The  domain  of  interest  is  many  times 
larger  in  the  principle  flow  direction  than  in  the  transverse  directions  which  indicates 


the  use  of  the  parabolized  form  of  the  Navier-Stokes  equations.  At  regions  very  close 
to  the  body,  reversed  flow  makes  these  techniques  invalid.  At  locations  fax  enough 
from  the  body  that  the  axial  velocity  is  unidirectional,  the  flow  field  can  be  adequately 
described  in  the  upstream  plane  only  by  using  all  velocity  components,  turbulent  model 
components  (such  as  k  and  c),  and  the  pressure.  For  most  numerical  schemes  this 
represents  an  over-specified  condition. 

Two  methods  are  shown  to  deal  with  this  problem  within  the  context  of  the  par¬ 
tially  parabolic  equations.  In  both,  the  primary  concept  is  to  separate  the  calculation 
used  for  local  mass  continuity  and  stream  wise  pressure  gradient.  In  both  cases  primary 
emphasis  is  placed  on  the  local  continuity  equation  and  the  results  of  that  calculation 
is  used  to  estimate  the  streamwise  pressure  variation. 

Of  the  two  methods,  the  procedure  which  uses  an  exact  equation  to  calculate 
streamwise  pressure  variations  and  an  approximate  transverse  continuity/pressure 
equation  appears  to  match  experimental  measurements  more  closely.  The  procedure 
which  uses  the  approximate  pressure  correction  method  for  both  streamwise  pressure 
variation  and  transverse  continuity  converges  much  more  rapidly,  requires  much  less 
computer  storage,  and  makes  predictions  reasonably  close  to  the  total  pressure  formu¬ 
lation  and  experiment. 

Both  methods  resut  in  good  quantitative  agreement  with  the  measured  velocity 
data  in  the  near-wake  region  and  do  so  without  the  necessity  to  empose  the  emperical 
pressure  data  as  was  done  in  the  QPP  calculation. 
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Table  1  —  Parabolic  Form  of  the  Governing  Equations 


Momentum: 


MU  x.  ^VU)  x.  d<WU')  -  ?<!±  x.  ?lL y  dP  x.  r 

P  dx  +  p  dy  +  P  dz  dy  +  dz  dx  +  F‘ 


d(UV )  .  3(V2)  .  d(WV)  Sr ‘  dr '  dp 

"ST  +  p-aT  +  p—z - dT  +  —  ~  + 


d(UW)  .  d(VW)  ,  d{W J)  3<„  dr'WLt  dP 

P~~5T'  +  P~W  p~dT  =  ~  ~sr  -  37  +  F- 

< 

Continuity: 


dU  dV  dW 
dx  ^  dy  ^  dz 


Turbulent  Kinetic  Energy: 


d(Uk)  d(Vk)  d(wk)  d  ,,ckk-j  .  dk  x  d  ,,ckk—  ak 

V  +  V  +  -  V  +  *«— ■*  -  ">3; : 1  +  ^  ~  *>a7l+ 


d  ,Ckk _ dk,  d  ,Ctk _ dk,  _ dU  _ dU 

— ( - V/W/—)  +  ■£-( - Vtwf—)  +  Uivt—  +  ufw/-- - h  e  =  0 

dy  e  dz  dz  e  dy  dy  dz 


Turbulent  Dissipation: 


%Ue)  ,  d(Ve)  ,  d(We)  ,  dfCkk-^de,  ,  d  tCkk  —  dt 

- dT  +  -dT  +  -dT  +  +  37(*7  u';  37 


d  .Ctk _ de .  d  Cek _ de  s  e_. _ <917 

— ( - Vtw(-)  +  — ( - Vtwf -T-)  4-  -C]utvt— 

dy  e  dz  dz  e  dy  k  dy 


+  jC]  u/wf  +  Cf  €—  =  0 

fc  ar  A; 


Where: 


r4  is  the  total  Reynolds  and  molecular  stress, e.g.  r4  r.  = 

is  calculated  according  to  which  method  is  used. 
P(y,z)  is  adjusted  to  satisfy  continuity 

u, x  act  in  the  axial  or  streamwise  direction 

v, y  act  in  the  horizontal,  transverse  coordinate  direction 

w, z  act  in  the  vertical,  transverse  coordinate  direction 
F  is  the  momentum  source  term 


{Ck,Ct,C},C*j  =  {.12,  .09, 1.44, 1.92} 


pu'v' 
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Table  1  (cont.) 


Reynolds  Stresses: 

■  +(%)*) -2c,£(g) 

W*  =  c,k  -  c2c,£(f  ?  -  2Ct$(%) 

^  =  c3k- c,c,£(%f  - 2C,£(eg) 

■  S7J7=-C^f -CJC.4j(t(^  +  ^)  +  2f(f +  1^)) 

Tww—  r  hLau  r  r  *3 1 au <  aw  ,  av  N  ,0au ,  au  ,  aw  ^ 
u/wt  -  -C4T  -57  -  C2C47T(a7(  97  +  a?)  +  2  37  (  37  +  -57)) 

77^7  _  ,  aw>v  rrki(dUdU\ 

u;/u/--C4T(a7  +  -ay-)-C2C4  7T(-57a7) 

{Ci,C2,C3,C4}  =  {.94,  .067,  .56,  .068} 


Pressure  Equations: 

Total  Pressure  Formulation: 

6 

=  PD'iP<A{  +  (p^u  ~  pUd)AyAz  +  (pFe  -  pKJA.*' Ar  +  {pWn  -  pWa  )A  xAy 
1 


where 

^  * 

U,  are  the  reduced  velocity  components 

c?.  =  3hE;^“+s»+fl,.."i 

^_=3hfE *,  Ai’ +  B' +  Dr,.’] 

w .  =  ^E)  A,-  +  B»  +  D,,.»] 

P,  are  the  pressures  at  location  i 

Ai,Bi,D{ ,  axe  coefficients  from  the  momentum  equations 

Pressure  Correction  Formulation: 

6 

A pP'  =  ^pD-P'Ai  +  (p!7u  -  pUd)AyAz  +  (pVe  -  pVw)AxAz  +  (p\Vn  -  p\V3) AxA/y 


where 

U,  are  the  calculated  velocities  at  the  control  volume  faces  before  continuity  correction. 
P[  are  the  pressure  corrections  at  location  i 
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Fig.  1  -  Axial  Distributions  of  Minimum  and  Maximum  Streamwise  Velocity 
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22.7%  MAXIMUM  SWIRL 


Fig.  3  -  Swirl  Vectors  at  x/D=  I 
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3.8%  MAXIMUM  SWIRL 
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.  4  -  Swirl  Vectors  at  x/D-28 
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Fig.  9  -  Axial  Distribution  of  Minimum  and  Maximum  Pressure 
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(c)  PRESSURE  CORRECTION  CALCULATION 
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